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ABSTRACT 

We investigate the migration of low-mass planets (5M 0 and 2OM 0 ) in accretion 
discs threaded with a magnetic field using 2D MHD code in polar coordinates. We ob¬ 
served that, in the case of a strong azimuthal magnetic field where the plasma parame¬ 
ter is (3 ~ 1—2, density waves at the magnetic resonances exert a positive torque on the 
planet and may slow down or reverse its migration. However, when the magnetic field 
is weaker (i.e., the plasma parameter [3 is relatively large), then non-axisymmetric 
density waves excited by the planet lead to growth of the radial component of the field 
and, subsequently, to development of the magneto-rotational instability, such that the 
disc becomes turbulent. Migration in a turbulent disc is stochastic, and the migration 
direction may change as such. To understand migration in a turbulent disc, both the 
interaction between a planet and individual turbulent cells, as well as the interaction 
between a planet and ordered density waves, have been investigated. 

Key words: accretion, accretion discs — magnetic fields — MHD — waves — planets 
and satellites: dynamical evolution and stability — planet-disc interactions 


1 INTRODUCTION 


The migration of protoplanets embedded in accretion discs 
has been extensively studied in recent years (see, e.g., review 
by jKley & Nelson|2 012 ). A low-mass protoplanet excites den¬ 
sity waves in the disc at the Lindblad resonances and migrates 
toward the star if no other torques are exerted (e.g., |Goldreich| 
|& Tremainej 19 79| ). However, a corotation torque originating 
near the planet can slow or even reverse the migration (e.g., 
jTanaka et al.|2002]|Masset|2002HMasset et al.|2 006). Another 
possible cause of slowed or reversed migration is the exci¬ 
tation of magnetic resonances in the disc when a relatively 
strong (/3 ~ l[j) azimuthal magnetic field is present (Terquem 


2003 ). These magnetic resonances can halt or reverse the mi 
gration of the planet if the magnetic field is sufficiently strong 
and it has a steep gradient toward the star. The action of this 
mechanism has been shown in numerical simulations by |Fro-| 
jmang et~al] ( [2005] ) . 


* E-mail of corresponding author: mcomins@astro. Cornell. edu 
1 (3 = 87t P / B 2 , where P and B are the pressure and magnetic field 
in the disc, respectively 


On the other hand, accretion discs threaded with a weak 
magnetic field are prone to the magneto-rotational instability 
(MRI; Balbus & Hawley 1991, 1998 ). Numerical simulations 
of planet migration in a turbulent disc have shown that the 
semimajor axis of the planet changes stochastically due to the 
planet’s interaction with the turbulent cells in the disc. The 
migration rate can increase or decrease, and the direction of 
the migration can reverse ( jNelson & Papaloizou[ 2004 ). 

An initial goal of this paper was to more deeply inves¬ 
tigate the effect of magnetic resonances on the migration of 
a low-mass embedded planet following the methods used by 
|Fromang et al.| ( |2005| , albeit that we explored a larger pa¬ 
rameter space with regards to the mass of the planet and 
the density distribution in the disc, as well as the magnetic 
field strength and its distribution in the disc. While exploring 
this expanded parameter space, we noticed the development 
of MRI-driven turbulence in the disc in many cases. In par¬ 
ticular, when the magnetic field is weak, non-axisymmetric 
density waves excited by the planet lead to perturbations in 
the initially-azimuthal magnetic field and, subsequently, to 
growth of the radial component of the magnetic field near 
the planet. The azimuthal field grows with time due to the 
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stretching of the radial field lines by the differential rotation 
of the Keplerian disc. As a result, an MRI-type instability de¬ 
velops and propagates to larger distances from the planet, and 
the disc becomes turbulent. The planet interacts with the tur¬ 
bulent cells, and the torque associated with this interaction is 
much larger than the torque seen prior to the development 
of turbulence (i.e., the torque due to magnetic resonances). 
Therefore, we also studied the parameter space describing 
the transition from a laminar to a turbulent disc, as well as 
the migration of a planet in a turbulent disc in detail. 

We first describe how the migration rate and direction of 
a low-mass planet are affected by the planet mass, as well as 
the distribution of surface density and magnetic field in the 
laminar disc. We also show how these parameters alter the 
time at which the onset of turbulence occurs in magnetized 
discs. Finally, to better understand the interaction between 
the planet and turbulent cells in the disc, we study the torque 
on the planet by (1) an ordered, low-amplitude density wave 
generated at the inner boundary and (2) a high-amplitude 
wave that is excited in a turbulent MHD disc. 

The plan for this paper is as follows. In ^2] we overview 
the theory of the different sources of torque on the planet. 
In ^3] we describe our physical model and numerical setup. 
In ^fjwe describe our main parameters. We describe our test 
simulations in the hydrodynamical disc in ^5] and simulations 
in a laminar magnetic disc in ^6] Migration in turbulent discs 
is shown in ^7] and migration under the influence of density 
waves is shown in 5j§] We conclude in ^9] 


2 THEORETICAL BACKGROUND 

In this section, we discuss the causes of inward and outward 
migration in hydrodynamic and MHD discs. 


2.1 Migration in hydrodynamic discs 


In the absence of a magnetic field in the disc, the migra¬ 
tion of the planet is determined by the Lindblad and coro¬ 
tation torques ( [Goldreich & Tremaine| 1 19791 |Ward| [1986] 
|1997| ). The Lindblad torque is generated when the planet 
excites ra-armed waves in the disc with “orbital” frequen¬ 
cies lj = mQ p . The dispersion relation for these waves in 
a low-temperature, low-mass disc is (u — mQ) 2 — k 2 , where 
k, = f^Kep is the epicyclic frequency in a Keplerian disc. Sub¬ 
stituting O(r) = y / GM*/r 3 , the locations of the Lindblad 
resonances are 


r*LR 




mil 

m 


2/3 


( 1 ) 


For each value of m (except m — 1), there is one resonance 
located closer to (farther from) the star than the planet called 
the inner (outer) Lindblad resonance. The inner Lindblad res¬ 
onances exert a positive torque on the planet and “push” the 
planet outward, while the outer Lindblad resonances have the 
opposite effect. The total Lindblad torque (also known as the 
“differential Lindblad torque”) determines the migration di¬ 
rection in the absence of other torques. The differential Lind¬ 
blad torque is negative, resulting in overall inward migration 
of the planet if there are no other sources of torque. 

However, there is typically a corotation torque exerted on 
the planet as well. The physics of the corotation torque and 


its effect on planet migration has been studied by number 


of authors (e.g., [Paardekooper & Mellema 2006; 

Baruteau & 

Masset 2008; 

Paardekooper & Papaloizou 2008; Kley et al. 

2009 

Masset & Casoli 2009 , 20lOj; Paardekooper 

et al. 2010 

2011 

). In an isothermal disc, both the magnitude and sign of 


the corotation torque depend on the gradient of the surface 
density at the corotation radius. 


2.2 Migration in magnetized discs 

2.2.1 Migration in laminar MHD discs due to magnetic 
resonances 


In strongly magnetized laminar discs, magnetic waves can be 
excited, and torques associated with the magnetic waves can 
affect a planet’s migration, (e.g., Terquem |2003[ [From ang 
|et al.|[2005{ |Fu & Lai||2011| ). |Terquem ( 2003| ) investigated 
the propagation of waves in a magnetized disc in which the 
magnetic field is purely azimuthal and found that there are 
two singular radii at which the frequency perturbation in a 
frame rotating with the fluid matches that of a slow MHD 
wave propagating along the field lines; these radii define the 
locations of the magnetic resonances. The inner and outer 
magnetic resonances are denoted similarly to Lindblad res¬ 
onances, with the inner magnetic resonances hmr < r p , and 
the outer magnetic resonances tomr > r p . Terquem ( 2003 ) 
derived the following dispersion relation: 


m 2 (Vl — Q p ) 2 — 


m 2 c 2 v a 


( 2 ) 


Here, va is the Alfven speed, given by |Terquemj (2003) as 
v\ — (B 2 )/4 ttE, where ( B 2 ) = / B 2 dz is the vertically- 
integrated square of the magnetic field. The sound speed 
c 2 = d(P)/dE, where ( P ) = / Pdz. |Terquem| ( 2003) de¬ 
termines the strength of the magnetic field via f3 Va =~<? s /v 2 A , 
which is evaluated at the location of the planet. [^] 

The locations of the resonances when the disc is Keple¬ 
rian (fl « Qk) and thin ( H/r <C 1) are given by 


| r M - r p \ = 


2 H 

3\/l + Pv A 


(3) 


where the thickness of the disc H and the plasma parameter 
/ 3 Va = (? s /v\ are evaluated at r — r p . As the field becomes 
weaker (i.e., p VA oo), the magnetic resonances converge 
toward the corotation radius. We suggest that c s /(rQ) « H/r 
and find the position of the inner and outer magnetic reso¬ 
nances to be 


n mr — r p 


2 H 

3 a/1 + (3 va 


T OMR = r p + 


2 H 

3a/1 + Pv A 


(4) 


jTerquem] ( 2003) showed that the waves associated with the 
magnetic resonances can exert a positive torque on the planet 
that is larger in magnitude than the differential Lindblad 
torque if the magnetic field increases toward the star. 

|Fromang et al.| ( |2005| ) performed simulations of a 
planet’s migration in a magnetized disc and found that the 
magnetic resonances can slow or stop the migration of a 
planet. In particular, they investigated the migration of a 5 M© 
planet embedded in a disc that is threaded with an azimuthal 


2 In our simulations, we use the standard definition of the plasma 
parameter, ft = 2 (3 VA . 
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magnetic field with B. f oc r~ k with E = const initially. They 
took an initial value of /3 Va = 2 at the planet’s location and 
found that, depending on the steepness of the magnetic field 
distribution (defined by k), the planet’s migration slowed or 
reversed in direction. 

|Guilet et akl ( |2013| ) studied the effects of an MHD coro¬ 
tation torque on a low-mass planet in a 2D laminar disc with 
a weak azimuthal field threading the disc. The field was not 
strong enough to generate an appreciable torque from mag¬ 
netic resonances, and it was not strong enough to dominate 
the hydrodynamic corotation torques, but a “torque excess” 
attributed to the presence of the magnetic field was found. 


small compared with r), oc exp (ik r r + imtp + ik z z — iut), un¬ 
der the conditions k 2 z » k 2 and ( k^VA ) 2 < -rdQ 2 /dr, where 
km m/r. This is the same as the condition for the MRI in¬ 
stability in a disc threaded by a vertical field with B z B^ 
and k z —>> k^. Therefore, perturbations of the azimuthal field 
may also lead to MRI turbulence. 

Work done by |Baruteau et al.| ( |2011| ) and |Uribe et ak] 
( 120111 ) showed the existence of additional MHD corotation 
torques in MRI-turbulent discs using 3D MHD simulations. 
Furthermore, Nelson & Papaloizou (2004) investigated planet 
migration in an MRI-turbulent disc and observed stochastic 
migration. However, the origin of the torque on the planet in 
a turbulent disc is not well studied. 


2.2.2 Migration in a turbulent disc 


The MRI arises under conditions in which a weak magnetic 
field threads the disc and the radial component of the field 
can be stretched and enhanced by the differential rotation in 
the disc. Such a seed field can either be an axial field, perpen¬ 
dicular to the disc ([Balbus & Hawley|1 99l]|l998 ), or an az¬ 
imuthal field (Terquem & Papaloizou 1996). Below, we briefly 
summarize the conditions for the onset of the MRI instability 
in these two cases. 

Axial Seed Magnetic Field. 

Balbus & Hawley ( |1991| ) considered the case in which an 
axial magnetic field, B 0 z, threads a Keplerian disc that ro¬ 
tates with an angular speed Q. For axisymmetric perturba¬ 
tions of the disc, for which 5v = [Sv r (z,t),Sv (p (z,t) 1 0] and 
£B = [SB r (z,t),SB (p (z,t), 0], and for perturbations propor¬ 
tional to exp (ik z z — iut), one finds the dispersion relation 


LU± ( k z VA ) 2 + 




+ 


4 (k z VA&) 



where va = Bo/y/^Tip is the Alfven velocity, and 
K r = [4fl 2 + 2 rQdQ/dr] 1 ^ 2 


(5) 

( 6 ) 


is the radial epicyclic frequency of the disc. In order for the 
perturbation to fit within the vertical extent of the disc, one 
needs k z h > 1, where h — c s /Q is the half-thickness of the 
disc and c s is the isothermal sound speed in the disc. For most 
conditions, the disc is thin, with h <C r or c s <C rQ. 

Evidently, instability can occur if uj 2 _ < 0, which happens 
if ( k z VA ) 2 < —rQdQ 2 /dr. For a Keplerian disc, this corre¬ 
sponds to ( k z VA ) 2 < 3f) 2 . Therefore, the above-mentioned 
condition that k z h> 1 implies that instability occurs only for 

VA < C s . 

Toroidal seed magnetic field. 

|Terqu em & Papaloizou |([l 996) studied the linear MHD stabil¬ 
ity/instability of a thin Keplerian disc with a toroidal magnetic 
field B^ (r, z). This case is more complicated than the case of 
a vertical field ( [Balbus & Hawleyfl 998 ). 

The complication in the toroidal field case results from 
the presence of both the MRI instability and the buoy¬ 
ancy instability of the toroidal field ( Hoyle & Ireland 1960; 
Parker 1966). The buoyancy instability is triggered by an 
azimuthally-dependent, radially-localized vertical displace¬ 
ment in the plasma and the toroidal field (Terquem & Pa-| 
paloizou 1996 ). These authors find that the MRI instability in 
a thin disc with an embedded toroidal magnetic field shows 
up in localized perturbations (i.e., whose wavelengths are 


3 MODEL 

3.1 MHD Equations 

We utilize the MHD equations to numerically evaluate the 
perturbative effect of the planet on the disc: 


(i) Continuity equation (conservation of mass) 

^ + IA(rE« r )+IA( Ew ) = 0 , (7) 

at r or r dtp 

where E = f pdz is the surface density (with p the volume 
density), and v r and are the radial and azimuthal veloci¬ 
ties, respectively. 

(ii) Radial equation of motion (conservation of momen- 


turn) 




+ 

1 d 

r dr 


+ 

1 d 


r dtp 


— 

n 

- + 


T.vl + n + — 


4/ rep 

ZJ V r V<n -— 

47r 




8 n 


4/ rr 

47T 


^rr + ^ 

87rr 


w ^ GM+ 

— - z]---b z] Wr 


( 8 ) 


where n = / Pdz is the surface pressure (with P the volume 
pressure), Ew r is the radial force exerted on the disc by the 
planet (per unit area of the disc), and Tf rr , 4^, and 4^ are 
magnetic surface variables that are discussed in more detail 
in$m] 

(iii) Azimuthal equation of motion (conservation of angu¬ 
lar momentum) 




}_d_ 

^ r 2 dr 
1 d_ 
r dtp 


EiV^Vip 


T/ 


r <p 


+ -7V- 


P'Vip + n + 


47T 

^rr + 

8tv 


T/ 


47T 


(9) 


where Ew^ is the azimuthal force exerted on the disc by the 
planet (per unit area of the disc). 

(iv) Radial induction equation 


d<F r 1 d 

-W + rdf^~ V ^ ) = °’ 


( 10 ) 


where 4> r and are magnetic surface variables, also dis¬ 
cussed in j ]3.1.1| 

(v) Azimuthal induction equation 


d$ v , d f _ 
~M + d~r 


V v $r) = 0. 


( 11 ) 
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(vi) Entropy balance equation 

I (ES)+ 11 {rESvr) +1 ip {ESv ^=°’ (12) 

where S = n/E 7 is a function analogous to entropy and 
we use 7 = 1.01 so that our disc is isothermal. We chose an 
isothermal disc to ease comparisons of our results with the 
results of other authors (e.g., by |Fromang et al.| 2005). 

Moreover, viscosity terms were added to the equations 
of motion following the a prescription of Shakura & Sun- 
|yaev| ( |1973| ) (see details in [Koldoba et al.||2015|). We use a 
very small viscosity (a = 0.001), analogous to|Fromang et 
[ah] ([2005 ), to isolate the effects of slow magnetosonic waves 
in the disc by smoothing out the effects of fast magnetosonic 


and the azimuthal equation of motion becomes 




+ 7A 

o o 


J2 l 


r YiVrVtn 


177 


1 d 7 2 „ + 

+ s< + n + 

r c)(f \ 


( p 


87 r 


$ 2 

_ jp 

47T 


( 18 ) 


The radial and azimuthal induction equations become, re¬ 
spectively, 


and 

~dt^ + dr ~ v^r) — 0 . 


(19) 


( 20 ) 


3.1.1 Magnetic surface variables 

The “volume” values for the radial and azimuthal magnetic 
fields are given by B r and B ip . Their vertically-integrated 
counterparts may be defined similarly to E and II as 


-/ 


4> r = / B r dz 


and 


‘I’,: = 


-/ 


B^dz, 


(13) 


respectively There are also terms in the MHD equations in¬ 
volving magnetic flux or energy that involve products of 
these variables: T? 2 , T? 2 , and B r Bf. We define the following 
vertically-integrated quantities for these products, 


4 >rr — J Brdz , 4^ — J B^dz, and 47^ = j B r B^dz, 

(14) 


respectively 

As shown in Eqns. ©-GD. the induction equations use 
4> r and 47 , while the equations of motion use ^ rr , 47 ^, and 
4^. As such, we need a way to relate <f> and 47 We can do 
this by introducing a “magnetic” thickness of the disc, hm. 
Using the definitions of <f> and 47 with hm, we find that 


47r 


4> r 4> r 

Hm 


^ hm ’ 


$ - 77 

h m ■ 


(15) 


We suggest that the value of hm = const is the same in all 
three relations. By relating 4> and 4^ in this way, we can de¬ 
fine a “surface magnetic field”, 03, such that the MHD equa¬ 
tions are parameterized with respect to a single magnetic field 
variable, 


03 r = 


4> r 

a/ hm 


and 


<B = -7_ 


(16) 


Then, the magnetic terms in Equations ([8]) - ([TlJ take their 
usual form. 

The radial equation of motion becomes 

>2 1 ^2 N 

8tv 47t 


s' 21 ") 


+ -77- 


1 d_ 
r dr 


vt + n + * + ** 


1 a / SrS, 

+ r 47T 


YiV\ 


- E 


2 

r 

GM , 


n 03 2 + 03 2 

r 87rr 

+ E 1/7 , 


3.2 Planetary equation of motion 

We calculate the equations of motion in the stellar reference 
frame, which is not inertial because the star also revolves 
about the center of mass of the system. So, an inertial force 
term is added to the equation of motion for both the disc and 
the planet. Assuming that the inertial acceleration is only due 
to the gravitational attraction between the star and the planet 
(but not the disc), the inertial force per unit mass (i.e., accel¬ 
eration) is 

GM P 


r 3 A P- 


( 21 ) 


To describe the gravitational influence of the planet on the 
disc, we use a gravitational potential similar to that used by 
|Fromang et al.| ( [2~005| ), 




GM V 


y/r 2 + r 2 _ 2 rr p cos((^ - tp p ) + < 


( 22 ) 


where e = 0.1 H is the gravitational smoothing length, and H 
is the scale height of the disc. The total force per unit mass is 


W = Wp+Wj, W P = -V4>p. 


(23) 


The components of this force in polar coordinates, w r and w 
are used in Eqs. and ( [T8] ). 

The force exerted on the planet by a particular fluid ele¬ 
ment with mass dM = T^rdrdtp, is the acceleration given in 
Eqn. ( [23] ), with opposite sign, multiplied by the mass of the 
fluid element, 


dfdii 


= —dMw v — dMV4 D . 


(24) 


We then calculate the total force exerted on the planet by 
the disc by integrating over the disc within the computational 
domain, 


^ disc—>p 


— / dfdisc—>-p — / 

J disc J dis 




(25) 


/disc «/disc 

which we use to find the position (r p ) and velocity (v p ) of the 
planet at each time step via the planet’s equation of motion: 




GM*M P 

73 

1 p 


GMi 


H - F 


(26) 


The gravitational torque in the 2 direction on the planet is the 
sum over the torques from each fluid element: 


(17) '-hz — / [r x dfdisc—>-p] 

J disc 


(27) 
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We also calculate the planet’s orbital energy and angular mo¬ 
mentum per unit mass (e.g., Murray & Dermott 1999) via 

■^=^I v p| 2 -——- and L = r p xv p , (28) 

Z Tp 

respectively. We use these relationships to calculate the semi¬ 
major axis and eccentricity of the planet’s orbit at each time 
step, 


a— — 


1 GM* 

2 E 


respectively. 


and 


e = 



L 2 

GM*a 


(29) 


3.3 Initial conditions 


The initial conditions for the surface density surface pressure 
and surface “magnetic field” are defined to be power laws 



(30) 

(31) 

(32) 


where r* is a characteristic radius in the disc; in most of our 
simulations n = r p ,z is the initial location of the planet. Addi¬ 
tionally n, l, and k are the power laws in the density pressure, 
and magnetic field distributions, respectively. We suggest that 
the initial pressure distribution is similar to the initial density 
distribution (i.e., n — l). The initial sound speed in the disc at 
n, defined as a fraction of the Keplerian speed at this radius, 
k s , is 


2 


C* 


GM* 

n 


(33) 


where k 8 — 0.01. The initial surface pressure at r* is 
Hi = EzC^. The initial plasma parameter at n is then de¬ 
fined as 


Pi = 


87r]Tj 
<B 2 ■ 


(34) 


Note again that this plasma parameter, ft, is twice as large as 
the plasma parameter defined in Froman g et al.| ([2005). We 
determine the surface magnetic field distribution via 


Variable Definitions 

ro 

Reference distance 


r 0 

Reference velocity 


Po 

Reference orbital period 



Reference disc mass 


E 0 

Reference surface density 

Bo 

Reference magnetic field strength 

To 

Reference torque per unit mass 

(a) 

Standard reference values 


r 0 = 0.1 AU 

r 0 = 1 AU 

r 0 

94.1 km s -1 

29.8 km s -1 

Po 

11.6 days 

367 days 

M d0 

1 Mq 

1M© 

E 0 

8.84 x 10 8 g cm -2 

8.84 x 10 6 g cm -2 

Bo 

228 kG 

2.28 kG 

To 

8.85 x 10 13 cm 2 s“ 2 

8.85 x 10 12 cm 2 s“ 2 

(b) 

Rescaled reference values 


r 0 = 0.1 AU 

ro = 1 AU 

ro 

94.1 km s -1 

29.8 km s —1 

Po 

11.6 days 

367 days 


io- 3 m© 

10- 3 M© 

So 

8.84 x 10 5 g cm -2 

8.84 x 10 3 g cm -2 

Bo 

7.22 kG 

72.2 G 

To 

8.85 x 10 13 cm 2 s“ 2 

8.85 x 10 12 cm 2 s“ 2 


(c) 


Table 1. Example reference units calculation. The stellar mass is 
M* = 1 Mq. (a) Definitions of the variables used, (b) Standard refer¬ 
ence values, corresponding to S = 1. (c) Rescaled reference values, 
corresponding to E = 0.001, which is the value used in the code. 


ft = 100 to study migration in a turbulent disc. We consider 
the migration of a 5M® planet in most of the simulations, as 
well as a 20M® planet in a few simulations. 


3.4 Reference Units 


93c® i — 



(35) 


We determine the initial equilibrium in the disc from the 
force balance in the radial direction: 


_< 1 

r E dr 87rr 2 E 


d(HB <ft 2 GM* 

' o 


dr 


= 0. 


(36) 


It is satisfied if the initial distribution of the azimuthal velocity 
has the form 


v 


2 _ 
<P ~ 



n 



1 — k 


Pi 



• (37) 


For all of the simulations presented here, we took 
k s — 0.01. The thickness of the disc at n is Hi = The 

smoothing radius of the gravitational potential is e = 0.1 H. 
An initial plasma parameter of ft = 1,2 is taken for analysis 
of migration in laminar MHD discs, while we increase ft up to 


Our simulations are performed using dimensionless units 
A = A/A 0 where A 0 are the reference units. We first choose 
some reference distance, ro. The results of our simulations 
are applicable to multiple regions in a protoplanetary disc, 
because ro can be chosen to correspond to different regions 
of the disc. In Table [l] we show examples of reference values 
for scale distances of ro = 0.1 AU and ro = 1 AU. The thick¬ 
ness of the disc is if = 0.1r o . We take the mass of the star, 
M * = 1 Mq, and determine the reference velocity, which is 
the Keplerian velocity at ro, vo = y / GMfto. The reference 
time is to = ro/vo, and we use the reference orbital period 
P 0 — 27rto as the unit of time in our plots. 

We next define the reference mass of the disc Mdo‘. 
Mao = M*. We then introduce the reference surface density, 
E 0 , such that M d o = E 0 ro. When the disc is homogeneous 
(i.e., E = const), Mao is the mass of the disc inside r = ro. 
The reference surface pressure is n 0 = E 0 ro. 
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The reference surface magnetic field is derived from the 
condition E 0 ^o = 53 q. Hence, 5$ 0 = y/^o^o- Taking into ac¬ 
count our definitions of the surface magnetic field (see Eqns. 
^3]and EH we obtain the reference volume magnetic field: 
Bo — %$o/y/ro- The reference torque per unit mass is defined 
as T 0 = r*o/t§. These reference units are used to convert equa¬ 
tions from dimensionless units. We show examples in Table 
®>. 


In our model, we use a small value of the dimensionless 
density E = 0.001; as a result, our dimensional characteristic 
disc mass is Mao = 10 _3 M 0 . The other characteristic values 
are closer to realistic values and are much smaller than the 
reference values shown in Table [T]d. We show in Table [TJ: the 
typical dimensional values corresponding to our simulations. 
One can see that, in the case of r*o = 1 AU, the surface density 
and other parameters are close to those expected in real pro- 
toplanetary discs. The values shown for r d = 0.1 AU are too 
large, but we keep these model parameters for all distances 
in order to compare our results with those in jFromang et al.| 
(2005) and others|^From this point forward, we will use only 
dimensionless units and remove tildes from variables. 


3.5 Grid and boundary conditions 

We use a polar grid that is uniform in the tp direction. In the 
r direction, however, the grid is non-uniform; the size of the 
grid cells increases with radius such that the sides are ap¬ 
proximately square-shaped throughout the disc. N r is defined 
to be the number of radial grid cells, and is the number 
of azimuthal grid cells. Our grid consists of N r — 480 and 
N(p = 1200 cells. Our inner boundary is at n n = 0.4r p ,i, and 
our outer boundary is at r out = 5 r Pji . 

We apply a wave damping procedure near both bound¬ 
aries similar to that used in §3.1.3 of Froman g et al.| ( [2b05| ) 
to avoid high-amplitude wave reflections that can over¬ 
whelm the migration signal. We perform damping after ev¬ 
ery time step for r < r damPjin , and r > r d am p,out, where 
r’dampjin — 1.375'Tm and 'TdampjOut — 0.8r*out* VVe calculate 
the velocities (y r , ry), as well as E, II, and S, using a method 
similar to Fromang et al. ( 2005), 


J * 


tTn T (»T »Tn) exp 
J/out T («T JTout) exp 


( r r damp,in \ 

l 5 in J 

{ r r damp,out \ 
y ^out J 


T ^ r’dampjin 


r* r’damp,out 
otherwise, 
(38) 


where J — (v r , v^, E, II, S), = 0.875 n n , and 

Sout = 0.8r out . Furthermore, Ji n and J ou t are the values 
of J at the inner and outer disc boundaries, respectively. 


3 It is often the case that the disc mass and the surface density are 
taken to be larger than in realistic discs. This leads to more rapid 
migration and thus shorter computational times (see, e.g., [Armitage| 
|& Rice|2006| . 


Type 

M p 

(M®) 

n 

k 

ft 

r p ,i 

Name 

Hydro 

5 

-1 

— 

— 

1 

H5n-lrl 

Hydro 

5 

-0.5 

- 

- 

1 

H5n-0.5rl 

Hydro 

5 

0 

- 

- 

1 

H5n0rl 

Hydro 

5 

1 

- 

- 

1 

H5nlrl 

Hydro 

20 

-1 

- 

- 

1 

H20n-lrl 

Hydro 

20 

-0.5 

- 

- 

1 

H20n-0.5rl 

Hydro 

20 

0 

- 

- 

1 

H20n0rl 

MHD 

5 

0 

0 

1 

1 

M5nOkO/31rl 

MHD 

5 

0 

0 

2 

1 

M5n0k0/32rl 

MHD 

5 

0 

1 

2 

1 

M5n0kl/32rl 

MHD 

5 

0 

2 

2 

1 

M5n0k2/32rl 

MHD 

5 

0 

0 

10 

1 

M5n0k0/310rl 

MHD 

5 

0 

0 

100 

1 

M5n0k0/3100rl 

MHD 

5 

1 

0 

2 

1 

M5nlk0/32rl 

MHD 

5 

1 

1 

1 

1 

M5nlkl/?lrl 

MHD 

5 

1 

1 

2 

1 

M5nlkl/32rl 

MHD 

5 

1 

1 

10 

1 

M5nlkl/310rl 

MHD 

5 

1 

2 

2 

1 

M5nlk2/32rl 

MHD 

20 

0 

0 

2 

1 

M20n0k0/32rl 

MHD 

20 

0 

1 

2 

1 

M20n0kl/32rl 

MHD 

20 

0 

2 

2 

1 

M20n0k2/32rl 

Waves 

5 

0 

— 


2 

W5n0k0r2 

Waves 

5 

0 

0 

100 

2 

W5n0k0/3100r2 


Table 2. Simulation names and their respective distinguishing vari¬ 
ables. See a more detailed description in ^4| 


4 PARAMETER SPACE 

We calculated a number of models at different initial param¬ 
eters, which are described below (see also Table |2j: 

Simulation type: The simulations are performed in a hy¬ 
drodynamic disc or an MHD disc, or to study the interaction 
between the planet and waves in the disc. 

Planet mass: We explored two different planet masses - 

5M® and 20M®. 

Surface density exponent, n: This defines the initial surface 
density distribution in the disc, according to E oc r~ n . 

Surface magnetic field exponent, k: This defines the initial 
surface magnetic field distribution in the disc, according to 

55 (xr~ k . 

Matter-to-magnetic pressure ratio, j3 % : This defines the 
value of fa at the initial location of the planet. 

Initial planet location, r p y. This defines the initial orbital 
radius of the planet. 

The last column of Table[2]shows the names of the simulations 
presented in subsequent sections. These names are referenced 
in the figure captions and related discussion. 
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Figure 2. Left Panel: The variation in semimajor axis versus time for a 5M® planet embedded in a hydrodynamic disc for different values of n 
(where E oc r~ n ). The models shown include: H5n-lrl, H5n-0.5rl, H5n0rl, H5nlrl. Right Panel: A similar plot, but for a 2 OM 0 planet. The 
models shown include: H20n-lrl, H20n-0.5rl, H20n0rl. 



the planet migrates towards the star. However, when the den¬ 
sity decreases towards the star, n — — 1, the corotation torque 
is larger than the Lindblad torque, and the planet migrates 
outwards. In the case of a more shallow positive density dis¬ 
tribution, n = —0.5, the Lindblad and corotation torques are 
almost equal, and the cumulative torque is small. 

The migration due to both the Lindblad and corotation 
torques has been studied in semi-analytical models by |Tanaka| 
|et al.| ( |2002p . In two dimensions, they obtained the torque on 
the planet from the disc: 


r total, 2D = —(1.160 + 2.828a) 


Mp Tpflp 

M* c s 


S(r p )r^Q: 


A 

p> 


(39) 


Figure 1. An example of a hydrodynamic simulation of migration 
corresponding to the model H5n0rl after 10 orbits of the planet. The 
color background shows the surface density distribution. The m = 
1, 2 Lindblad resonances are indicated by the solid black circles. 


5 MIGRATION IN A HYDRODYNAMIC DISC 


where a defines the slope of the density distribution 
(E oc r -a ), and c s is the isothermal sound speed in the 
disc at the orbital radius of the planet. This implies that the 
total torque on the planet in a two-dimensional disc is zero 
when a = —0.41, suggesting that 


As a first step, we investigated the migration of a planet in 
a hydrodynamic disc with an initially homogeneous density 
distribution (n = 0); these simulations are used as a base for 
studying migration in magnetic discs. Fig.[l]shows the surface 
density distribution in the disc after 10 orbits of the planet. 
The planet excites two density waves at the inner and outer 
Lindblad resonances (which are shown as black circles in the 
figure). 

As a next step, we considered different initial density dis¬ 
tributions in the disc, E ~ r~ n including those where the 
density is homogeneous, n — 0, increases towards the star, 
n— 1, or decreases toward the star, n = -0.5, -l.Fig.[2] (left 
panel) shows the variation of the planet’s semimajor axis over 
time for these cases. We observed inward migration when 
n — 0 and more rapid inward migration when n = 1. When 
n = —0.5, the migration almost stalls (where only slow out¬ 
ward migration has been observed), and we observed more 
rapid outward migration when n — — 1. 

The migration rate and its direction are determined by 
the cumulative value of the Lindblad and corotation torques, 
as described in Sec. |2.1| When the density increases towards 
the star (n = 1) or the density is constant in the disc (n = 0), 
the Lindblad torque is larger than the corotation torque, and 


(i) the torque on the planet is positive, and the planet mi¬ 
grates outward, when a < —0.41; 

(ii) the torque on the planet is zero, and the planet’s mi¬ 
gration halts, when a = —0.41; and 

(iii) the torque on the planet is negative, and the planet 
migrates inward, when a > —0.41. 


The value of n = —0.5 obtained from our simulations is very 
close to the value n — —0.41 corresponding to zero torque 
derived by Tanaka et al. (2002) for two dimensions. 

We also performed simulations of the migration of a 
2OM 0 planet. Fig. [ 2 ] (right panel) shows that the planet mi¬ 
grates inward when n — 0; almost no migration is observed 
when n — — 0.5 (as with the 5M 0 planet); and the migration 
direction is outward when n = — 1. These results are in accord 
with that of a 5M 0 planet, while the migration rate is faster 
in the 2OM 0 case. This is expected according to Eqn. ([39]), 
which states that the total Lindblad torque increases with the 
mass of the planet. 

Our simulations of hydrodynamic discs are in accord with 
theoretical studies (e.g., [Tanaka et al.[2 002 ) and simulations 
performed by others (see review by Kley & Nel son|2012| ). 
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Figure 3. The surface density and surface azimuthal magnetic field distributions for the model M5n0k0/32rl at t = 10. Left Panel: The two- 
dimensional surface density distribution. The dashed red lines show the positions of the magnetic resonances. Middle Panel: The two-dimensional 
surface azimuthal magnetic field distribution. The dashed black lines show the positions of the magnetic resonances. Right Top Panel: The one¬ 
dimensional surface azimuthal magnetic field distribution, taken along the horizontal dashed line in the left two panels. 

Right Bottom Panel: The one-dimensional surface density distribution, taken along the horizontal dashed line in the left two panels. The vertical 

dashed line shows the location of the planet 


6 MIGRATION IN LAMINAR DISC DUE TO MAGNETIC 
RESONANCES 

6.1 Magnetic resonances in the case of a constant 
density distribution 

We performed a series of simulations to investigate the influ¬ 
ence of an ordered azimuthal magnetic field in the disc on a 
planet’s migration. We started by investigating migration in 
discs with a homogeneous surface density distribution (i.e., 
n = 0 ), with initial and boundary conditions very similar to 
those used by Fromang et al7| ( |2005| ) . In particular, we consid¬ 
ered an isothermal disc (i.e., 7 = 1 . 01 ) and set the strength 
of the magnetic field near the planet such that fdi — 2. Our 
goal was to see whether our code can reproduce the results 
presented by Fromang et al. (2005}. 

First, we investigated a homogeneous magnetic field 
(i.e., k — 0). Fig. [3] shows an example simulation after 10 
orbital periods of the planet. We observed that the planet ex¬ 
cites ring-like waves at which the azimuthal magnetic field is 
stronger and the surface density is lower than in other nearby 
regions in the disc. These locations correspond to the mag¬ 
netic resonances that are excited by slow magnetosonic waves 
propagating along the field lines. Fig. [3] shows two rings of 
lower density (left panel), corresponding to two rings of en¬ 
hanced azimuthal field (middle panel). The position of the 
magnetic resonances is similar to that found in the simula¬ 
tions by |Fromang et al.] ( |2005| ), and they correspond to the 
theoretical resonance locations predicted by |Terquem| ( |2003| ) 
(see the dashed-line circles in Fig. [3]). The right panel shows 
the linear distribution of B^ and E in the radial direction. 

According to the theory presented by |Terquem| ([2003), 
these waves exert a torque on the planet that can reverse the 
planetary migration direction if the magnetic field distribu¬ 
tion is steep enough. We varied the steepness of the magnetic 
field, k (see Eq. [32]), taking k — 0,1 and 2, and we obtained 
different migration rates. Fig. [ 4 ] (left panel) shows that, for a 
constant magnetic field distribution {k — 0), the planet mi¬ 
grates inwards but more slowly than in the purely hydrody¬ 
namic case. When k — 1, the planet slowly migrates outward, 


while the planet migrates outward more rapidly when k — 2. 
These simulations confirm the result presented by ( Fromang 
|et al.|2005} : in discs with a relatively strong azimuthal mag¬ 
netic field, the magnetic resonances can slow or reverse a 
planet’s migration. 


6.2 Migration due to magnetic resonances for different 
density distributions 

Next, we investigated the action of the magnetic resonances 
at different surface density distributions. First, we took a den¬ 
sity distribution such that the density in the disc increases 
towards the star, n — 1, and repeated the above simulations 
at k = 0,1,2. Fig. |4] (middle panel) shows that the planet mi¬ 
grates inward in all three cases and, therefore, that the torque 
from the magnetic resonances is small compared with the dif¬ 
ferential Lindblad torque. The migration rate in the case of 
k — 1 almost exactly coincides with the purely hydrodynamic 
case. For a very steep magnetic field distribution, k — 2, the 
accretion rate is only slightly slower than for/c = 1. Over¬ 
all, we conclude that, in the case of this density distribution 
{n — 1), the positive torque associated with the magnetic res¬ 
onances is not strong enough to overcome the negative differ¬ 
ential Lindblad torque. 

In another example, we considered a disc whose density 
decreases towards the star, n — —0.5. This situation is pos¬ 
sible, for example, at the inner edge of the disc where the 
expanding magnetosphere or erosion of the disc may push 
the inner disc away from the star (e.g., [Lovelace et al.|2 008). 
When n — — 0.5, the migration rate in the hydrodynamic case 
is very low (see Fig. [3] top panel) because the negative differ¬ 
ential Lindblad torque is approximately compensated by the 
positive corotation torque. For this surface density distribu¬ 
tion {n— —0.5), the positive torque associated with the mag¬ 
netic resonances leads to outward migration of the planet, 
even for a flat magnetic field distribution, k = 0 (see Fig. [ 4 ] 
right panel). For a steeper magnetic field distribution, k m 1, 
the planet migrates outward even more rapidly. It appears, 
then, that magnetic torques may play a significant role at the 
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Figure 4. Change in semimajor axis versus time for a 5M® planet embedded in an MHD disc for different values of k (where 03 oc r ~ k ) 
and ^ = 2. The result for a hydrodynamic disc is shown for reference. Left Panel: Result shown for n — 0, corresponding to models 
H5n0rl, M5n0k0/32rl, M5n0kl/32rl, and M5n0k2/32rl. Middle Panel: Result shown for n = 1, corresponding to models H5nlrl, M5nlk0/32rl, 
M5nlkl/32rl, M5nlk2/32rl. Right Panel: Result shown for n = —0.5, corresponding to models H5n-0.5rl and M5n-0.5k0/32rl. 



Figure 5. Left Panel: Variation of the semimajor axis of a 2OM0 planet under the influence of Lindblad and magnetic resonances for various 
values of initial surface magnetic field exponent, k, for the case n = 0 and Pi = 2, corresponding to models M20n0k0/32rl, M20n0kl/32rl, 
M20n0k2/32rl, and the hydrodynamic case H5n0rl. Right Panel: Same as in the left panel, but over a longer interval of time. 


disc-cavity boundaries, as well as other regions where the sur¬ 
face density in the disc is flat or decreases towards the star. 

We conclude that the rate and direction of migration are 
determined by the steepnesses of both the magnetic field dis¬ 
tribution and the surface density distribution. If the surface 
density increases radially towards the star (e.g., n = 1), then 
the magnetic resonances do not exert enough torque to drive 
outward migration for any value of k. By contrast, when the 
surface density in the disc decreases toward the star (e.g., 
n = —0.5), the torque from the magnetic resonances is large 
enough to drive outward migration for different values of k. 
When the surface density is constant (i.e., n — 0), the mag¬ 
netic resonances drive outward migration when k — 1,2, and 
inward migration when k = 0. 

6.3 Migration of a 2OM 0 planet 

The above analysis shows that the effect of an azimuthal mag¬ 
netic field on a planet’s migration strongly depends on both 
the surface density and magnetic field distributions in the 
disc. It is also interesting to investigate whether the magnetic 
resonances also depend on the mass of the planet. To investi¬ 
gate this issue, we simulated the migration of a 2OM 0 planet 
in a magnetized disc with a flat surface density distribution 
(n = 0) and different steepnesses in the distribution of the 
magnetic field {k = 0,1,2). Fig. [ 5 ] (left panel) shows that, for 
a constant magnetic field in the disc, k — 0, the migration 
rate is slower than in the hydrodynamic case. For a steeper 
magnetic field distribution, km 1, the inward migration is 
even slower, and the direction of the migration reverses at 
an even steeper field distribution, k — 2. These results are 
in general agreement with those obtained for a smaller-mass 


planet: the magnetic resonances are strong enough to reverse 
the migration. Note that, at k = 2, the outward migration of 
the more massive planet is much faster than the migration of 
the lower-mass planet. 

However, longer simulation runs have shown that, for a 
2OM 0 planet, the laminar stage of the disc does not last long, 
< 10 — 40 planetary orbits. At later times, the disc becomes 
turbulent. Fig. [ 5 ] (right panel) shows that, after a relatively 
brief stage of slow migration in the laminar disc, the migra¬ 
tion becomes stochastic in the turbulent disc. We also see such 
a transition to stochastic migration for a 5M 0 planet, though 
at much later times (£ tU rb ~ 80 when n = 0, k = 0, and 
Pi = 2). In both cases, the disc becomes turbulent, and the 
semimajor axis varies in time stochastically due to the planet’s 
interaction with turbulent cells in the disc. 

We found that the observed turbulence is an interest¬ 
ing phenomenon and investigated the migration of planets 
in turbulent discs. Earlier, such migration was studied by |Nel-| 
|son & Pap aloizou (2004) in 3D simulations; they observed 
that the migration becomes stochastic and the migration rate 
may be strongly modified or reversed due to interaction with 
turbulent cells in the disc. In this paper, we performed 2D 
simulations in polar coordinates, and we see similar stochas¬ 
tic migration of the planet due to interaction with turbulent 
cells in the disc. While 2D simulations of the MRI are some¬ 
what restricted, they allow us an opportunity to investigate 
the details of the interaction between the planet and inho¬ 
mogeneities in the disc. Below, we investigate migration in 
turbulent discs. 
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Figure 6. Left Panel: Variation of the semimajor axis of a 5M® planet for in models with different values of fa when n = 0 and k = 0, 
corresponding to models H5n0rl, M5nOkO/31rl, M5n0k0/32rl, M5n0k0/310rl, and M5n0k0/3100rl. .Right Panel: Same as in the left panel, but 
over a longer interval of time. 



Figure 7. Surface density distribution and magnetic field lines corresponding to model M5n0k0/3100rl at t = 10,19, 22. This figure demonstrates 
how initially-azimuthal field lines (left panel) acquire a radial component from non-axisymmetric matter flow near the planet (density waves) 
and start forming loops (middle panel). Later on, the process spreads to larger distances from the planet, forming an inhomogeneous distribution 
of matter and magnetic field in the disc (right panel). 


7 MIGRATION IN A TURBULENT DISC 

In this section, we investigate the transition from the laminar 
to the turbulent disc, the formation of turbulent discs, and 
the migration of a low-mass 5M® planet in discs with differ¬ 
ent strengths of the magnetic field. According to the theory of 
the MRI (see Sec. |2.2.2) , instability is expected in magnetized 
discs with fa > 1; that is, in discs in which the Alfven veloc¬ 
ity v a — B/faAnp is smaller than the sound speed c s (e.g., 
jBalbus & Hawley| 1 991, 1998; [Terquem & Papaloizou] l996 ). 

Below, we show the results of simulations at different initial 
values of fa (see Sec. |7. 1|) a nd the details of migration in the 
turbulent disc (see Sec. |7.2| ). 

7.1 Migration in a turbulent disc with different fa 

We studied the migration of the planet in discs with differ¬ 
ent initial plasma parameters, fa = 1,2,10,100. Fig. [6] (left 
panel) shows that the disc is initially laminar, and the planet 
migrates smoothly. This period of smooth migration is longest 
when fa = 1 and 2, but it is shorter for fa = 10 and even 
shorter at fa — 100. This is understandable, because the mag¬ 
netic field is more easily tangled by the disc matter when the 
matter strongly dominates, such as when fa — 10 and 100. 
The right panel of the same figure shows the semimajor axis 
of the planet at later times. One can see that the disc is still 
laminar when fa — l. In all other cases, the disc becomes tur¬ 
bulent. Note that the migration often changes direction from 


inward to outward and vice versa; this reflects the stochastic 
nature of the migration process in a turbulent disc, as dis¬ 
cussed below in Sec. l7.2| 


7.2 Migration in a turbulent disc with fa = 100 

In this section, we consider the development of the MRI, as 
well as the migration of a planet in a turbulent disc, in greater 
detail. As a base, we consider a model where the initial value 
of plasma parameter is large, fa = 100 (model M5n0k0/3100), 
so that the MRI instability starts easily. 

We observed that the origin of the radial component of 
the field, which is required for the instability, is in the fact 
that the planet excites non-axisymmetric waves in the disc. 
These waves lead to non-axisymmetric motion of the gas in 
the disc and, subsequently, to the formation of a radial com¬ 
ponent of the magnetic field, which is further stretched by the 
differential rotation in the disc. Fig. [7] shows how parts of the 
initially azimuthal magnetic field near the planet (left panel) 
acquire a radial component and are subsequently stretched 
by the differential rotation in the disc, eventually forming a 
loop (middle panel). Later, this process occurs at larger dis¬ 
tances from the planet, many more field lines are stretched in 
the radial direction (such that different inhomogeneities and 
loops form), and the disc becomes globally inhomogeneous 
and turbulent. 

Fig. i (left panel) shows that the disc consists of 
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Figure 8. Planet migration in an MRI-turbulent disc in the model M5n0k0/3100rl at t = 50. Left Panel: The surface density variation in the inner 
part the disc. The planet’s location and the m = 1 and m — 2 Lindblad resonances are shown via the black circles. Middle Panel: The variation 
of the magnetic field magnitude (£?tot) in the inner part of the disc. The planet’s location and the m = 1 and m = 2 Lindblad resonances are 
shown via the white circles. Right Panel: The variation of the plasma parameter (/%) in the inner part of the disc. The planet’s location and the 
m = 1 and m = 2 Lindblad resonances are shown via the white circles. 



Figure 9. Planet migration in the model M5n0k0/3100rl shown for t = 5,40, 60. Top Panels: The variation of the surface density (E) in the inner 
part of the simulation region. The location of the planet, and the m = 1 and m = 2 Lindblad resonances, are shown via the black circles. Second 
Row of Panels from the Top: The one-dimensional variation of the surface density along the red line in the panel above each respective plot. The 
vertical dashed black line shows the location of the planet. Third Row of Panels from the Top: The torques acting on the planet from the disc. The 
green line shows the torque on the planet from the disc where r < r p , the red line shows the torque on the planet from the disc where r > r p , 
and the blue line shows the total torque on the planet. The dashed black line marks zero torque. Bottom Panel: The change in the semimajor axis 
of the planet over time. 
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azimuthally-stretched turbulent cells. The middle panel of the 
same figure shows that the total magnetic field, B to t, becomes 
strongly inhomogeneous, with some regions having the origi¬ 
nal polarity and others reversing polarity The distribution of 
the plasma parameter, , shows that the simulation region 
splits into regions that are either magnetically- or matter- 
dominated (right panel). Therefore, in the MRI regime, the 
disc becomes strongly inhomogeneous both in the density and 
in the magnetic field. 

We investigated the density distribution, torque, and 
semimajor axis evolution of the planet in this case in greater 
detail. Fig. [9] (top panels) shows the surface density distri¬ 
bution in the disc at t — 5,40 and 60. At t — 5, the disc 
is still laminar and two Lindblad density waves are clearly 
seen. However, after t « 10, non-axisymmetric motions start 
to “tangle” the field lines of the weak magnetic field, and an 
MRI-type turbulence gradually develops. MRI-type turbulence 
is observed during the time interval 10 < t < 60. 

The middle column (t = 40) shows small-scale turbulent 
cells that persist for many orbits. However, later (at t = 60), 
the turbulent cells become larger in size because “islands” of 
stronger magnetic field also become larger in size, and often 
one or two main density waves form in the inner parts of the 
disc. The beginning of this process is seen in the right column 
of Fig. [9] We suggest that the finite life of the MRI turbulence 
and formation of these larger-scale waves may be connected 
with the 2D nature of our MRI turbulence. However, the low- 
amplitude MRI turbulence proceeds over long periods of time, 
which is sufficient to study the migration of a planet in the 
turbulent disc. Formation of large-scale waves is possible in 
realistic discs; we use these waves to study the interaction of 
a planet with waves in the disc in Sec. [8] 

The second row in Fig. [ 9 ] shows the ID density distribu¬ 
tion along the line connecting the planet and the center of 
the star at the same times as the top row. When the disc is 
laminar, we see only a small bump in density associated with 
matter accumulation near the planet. However, as the disc be¬ 
comes more and more turbulent, larger and larger variations 
in the surface density of the disc are observed. 

The 3rd row in Fig. [ 9 ] shows the torques acting on the 
planet at the same times as the top two rows. When the disc 
is laminar (left column), the planet migrates due to the exci¬ 
tation of density waves at the Lindblad resonances: the inner 
torque is positive and smaller than the negative outer torque, 
so that the total torque is negative (the blue line in Fig. [ 9 ]). 
In this case, the planet slowly migrates inward (see the slow 
variation of the semimajor axis up to time t ~ 30 in the bot¬ 
tom row). However, when the disc is turbulent (t = 40,60), 
we observe that both the inner and outer torques can be either 
positive or negative; they both vary rapidly and the resulting 
torque also changes sign. 

The magnitudes of these turbulent torques are much 
larger than those seen in the laminar case (note the torque 
magnitudes shown on the y -axes in the third row of Fig. [9j. 
The torques become stochastic and correspond to the inter¬ 
action of the planet with individual turbulent cells. The vari¬ 
ation of the sign of the net torque shows also that the total 
averaged torque may be either negative or positive (i.e., that 
the direction of migration may change). In the shown exam¬ 
ple, the average total torque is negative and the planet mi¬ 
grates inward overall. Note that, in other cases, the direction 


of the migration may be either inward or outward (see, e.g., 
Figs. [5] an d@- 


8 INTERACTION BETWEEN A PLANET AND WAVES IN 
THE DISC 

The interaction between a planet and turbulent cells in a 
disc is a complex process. The planet interacts with a set 
of turbulent cells gravitationally, but the Lindblad density 
waves are not homogeneous. The closest cell may contribute 
to the torque more strongly than more remote cells, and 
the torque becomes more stochastic. Additionally, the planet 
passes through individual cells, and each cell may exert a 
corotation torque on the planet. It is difficult to track the in¬ 
teraction of a planet with individual turbulent cells. That is 
why we developed conditions in which a planet interacts with 
ordered inhomogeneities (in the form of waves in the disc). 
We consider two types of waves: (1) low-amplitude ordered 
waves generated in a hydrodynamic disc by a force at the in¬ 
ner boundary (Sec. |8.1| ), and (2) high-amplitude waves that 
form at later times in simulations using an MHD disc (Sec. 

HU) 

8.1 Interaction between a planet and low-amplitude 
waves in a hydrodynamic disc 

To better understand the interaction between a planet and 
individual turbulent cells, we created a model in which an 
ordered density wave is generated at the inner disc bound¬ 
ary by a periodic force that decreases with the radius as r -3 . 
This force generates density waves with a small amplitude, 
for which the density contrast between the wave and the disc 
is small, about 5 — 7%. The density wave then propagates 
through the simulation region. We placed a planet at an ini¬ 
tial radius of r p ,i = 2, away from the action of this force at 
the boundary. 

We observed that the planet moves faster than the wave 
and it interacts differently with different parts of the wave. 
Fig. [To] (top panels) shows slices of the surface density at 
t — 12,13,14,15,16. We show the moments when the planet 
is located at the inner or outer edge of a wave (i.e., where 
the surface density either decreases or increases toward the 
star, respectively). The corotation torque is larger than the dif¬ 
ferential Lindblad torque when the density slope in the disc 
corresponds to that in Equation ([39]); that is, the slope of the 
density distribution is either positive or only slightly negative 
toward the star (see Sec. [5j . This is expected when the planet 
is at the inner edge of a wave. 

The middle row of Fig. [To] shows the torques acting on 
the planet; the times from the top panels marked with verti¬ 
cal dashed lines. The total torque has maxima at t — 12 and 
t — 15. The top panels show that, at these moments of time, 
the planet is located at the inner edge of the wave and, there¬ 
fore, that the positive torque acting on the planet is the coro¬ 
tation torque (which appears due to the positive density gra¬ 
dient at the inner edge of the wave). The bottom panel of Fig. 
[To] shows that the planet migrates inward overall, because 
the differential Lindblad torque dominates overall. However, 
at t = 12 and t — 15, the planet migrates outward due to the 
temporarily dominant corotation torque. 

At t m 13 and t — 16, the planet is located at the outer 
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Figure 10. Migration of a 5M® planet in a hydrodynamic disc in the model W5n0r2, for which low-amplitude ordered density waves are 
propagated through the disc, at t = 12,13,14,15,16. The bottom two panels mark these times with vertical black dashed lines. Top Panels: The 
surface density (E) variation in the disc; the small black circle shows the location of the planet. Middle Panel: The torque on the planet from the 
disc. The green line shows the torque on the planet from the disc where r < r p ; the red line shows the torque on the planet from the regions of 
the disc where r > r p , and the blue line shows the total torque on the planet. The horizontal dashed black line shows zero torque. Bottom Panel: 
The change in the semimajor axis of the planet over time. 


edge of the density wave, where the density increases towards 
the star. At these moments of time, the total torque is negative 
and the planet migrates inward. We suggest that, at the outer 
edge of a wave, a planet excites density waves at the Lindblad 
resonances, and the differential Lindblad torque drives overall 
inward migration. However, when the planet is at the inner 
edge of the wave, the positive corotation torque is large, and 
the planet migrates outward. 


8.2 Interaction between a planet and high-amplitude 
waves in an MHD disc 

In this subsection, we analyze the passage of a planet through 
waves with much higher amplitudes. These high-amplitude 
waves often form at later times in simulations of MRI- 
turbulent discs, where the non-axisymmetry of the gravita¬ 
tional potential leads to the formation of the inner density 
waves. We chose the model where the initial plasma param¬ 
eter ^ ss 100 and a planet is placed at r p ^ — 2. At later 
moments in time, after the wave forms, the planet migrates 
to a radius of r p w 0.9 — 1. Fig. 11 (top panels) shows the 
density distribution in the wave and the position of the planet 
for several representative moments in time, where the planet 
is located at the inner edge of a wave (left and right pan¬ 
els, t m 180.25,184.75), the planet is in the middle of a 
wave (t = 181.5), and the planet is in the low-density region 
(t — 183). The density contrast between the wave and the rest 


of the disc is ~ 70%, which is about 10 times higher than the 
low-amplitude waves discussed in the previous subsection. 

The middle and bottom rows of Fig.[lT]show the torques 
acting on the planet and the variation of the planet’s semi¬ 
major axis in time. The dashed vertical lines show the four 
moments in time corresponding to the top four panels. One 
can see that, at t — 180.25, the planet is located at the in¬ 
ner edge of the wave, where the density decreases towards 
the star, and the positive corotation torque is expected to be 
larger than the differential Lindblad torque. Indeed, the mid¬ 
dle row shows that the total torque is positive at this moment 
of time, and the planet migrates outward. 

At the second moment in time, t = 181.5, the planet 
is located in the middle of the density wave. The middle 
row shows that, at this moment, a strongly negative torque 
dominates, and the planet migrates inward. This large nega¬ 
tive torque is primarily comprised of the differential Lindblad 
torque. This torque is proportional to the surface density in 
the disc, and therefore it is large when the planet is located 
inside the high-density wave. The positive corotation torque 
is relatively small. At t — 183, the planet is located away from 
the wave in the low-density part of the disc, and both torques 
are small. At the last considered moment (t — 184.75), the 
planet is again at the inner edge of the wave, where the coro¬ 
tation torque is positive and larger than the differential Lind¬ 
blad torque, and the planet migrates outward. 

In this example, the density wave is an analog of a large 
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Figure 11. Planet migration in high-amplitude waves in the model W5n0k0/3100r2 for t = 180.25,181.5,183,184.75. In the two bottom panels, 
these times are shown by vertical black dashed lines. Top Panels: Surface density variation in the inner part of the simulation region. The large 
black circles show position of the Lindblad resonances, while the small black circle shows the position of the planet. Middle Panel: The torque on 
the planet from the disc. The green line shows the torque on the planet from the part of the disc where r < r p ; the red line shows the torque on 
the planet from the disc where r > r p . The blue line shows the total torque on the planet. The horizontal dashed black line shows zero torque. 
Bottom Panel: The variation of the semimajor axis of the planet over time. 


turbulent cell, where the total torque is either positive or neg¬ 
ative and acts onto the planet depending on the position of 
the planet relative to the wave. We expect that the interaction 
with smaller-sized turbulent cells is similar to the observed 
interaction with MHD waves, but that the duration of the in¬ 
teraction is shorter. As a result of such interactions, a torque is 
exerted on the planet during short intervals of time, and the 
planet’s semimajor axis varies stochastically under the action 
of these torques. 

Based on our observations with both the low- and high- 
amplitude waves we can schematically describe the interac¬ 
tion between a planet and a turbulent cell (see Fig.p~2]). When 
the planet is at the outer edge of the cell (the side that is far¬ 
ther from the star), the density is increasing toward the planet 
(n > 0), and so we expect the torque to be negative and the 
planet to migrate inward. Conversely, when the planet is on 
the inner edge of the cell (the side that is closer to the star), 
the density decreases toward the star (n < 0), and so we 
expect the torque to be positive and the planet to migrate 
outward. This inward-outward migration does not “trap” the 
planet, because the cell is moving away from the star, so this 
is a transient process. 

9 CONCLUSIONS 

We investigated the migration of a low-mass planet (5M©) 
in magnetized discs using a Godunov-type HLLD MHD code 


^ -To star -To star 



I increases toward star (n > 0) I decreases toward star (n < 0) 

Torque < 0 Torque > 0 

Inward migration Outward migration 

Figure 12. Schematic showing the interaction between a planet and 
a turbulent cell. 

in polar coordinates ( [Koldoba et al.|2 015). The initial surface 
magnetic field is azimuthal, with different radial distributions 
B^ rsj r ~ k , and different strengths determined by the initial 
plasma parameter at the planet’s location, /3i, which varied 
between fa = 1 and 100. We also varied the initial radial 
surface density distribution in the disc, E ~ r~ n . Our main 
conclusions are as follows: 

(i) In strongly-magnetized discs (/% = 1,2), and where the 
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density distribution in the disc is flat (n = 0), the planet’s mi¬ 
gration is strongly influenced by magnetic resonances, which 
are excited near the planet and exert a positive torque on the 
planet. The migration slows down when the magnetic field 
distribution is flat (k = 0); it slows down more strongly when 
the magnetic field increases towards the star (k = 1); and the 
migration reverses when the field steeply increases towards 
the star (k = 2). These results are in accord with theoretical 
predictions by Terquem] (|2003|) and simulations by Fromang 
|etall ( [2005l ). 

(ii) Compared with Froman g et al.| ( [2005] ), we investigated 
the effect of magnetic resonances on the migration of a planet 
in discs with different density distributions. We observed that 
the steepness of the density distribution strongly influences 
the rate and direction of the migration. When the density in¬ 
creases towards the star (n = 1) the planet migrates inward 
at any steepness of the magnetic field (k = 0,1,2) and the 
effect of the magnetic resonances is negligibly small. This is 
because, at larger density steepness, the negative differential 
Lindblad torque is much larger than the positive corotation 
or magnetic torques. In the opposite situation, when the den¬ 
sity in the disc decreases toward the star (n = —0.5), and the 
positive corotation torque almost balances the negative Lind¬ 
blad torque, the role of the positive magnetic torque becomes 
very significant. The planet migrates outward due to mag¬ 
netic resonances at all values of the magnetic field steepness 
(k = 0,1,2). 

Ciii) Experiments with a larger mass planet, 20M®, show 
that the action of the magnetic resonances is similar to that in 
the case of a lower-mass planet: the positive torque increases 
with the steepness of the field, as predicted by the theory. 
We also observed that, at k = 2, the more massive planet 
migrates outward more rapidly than the lower-mass planet. 
This is in accord with Eqn. (65) of |Terquem] (2003]), in which 
the migration time scale is shown to be oc Mp. 

(iv) We investigated weakly-magnetized discs with initial 
plasma parameter fa = 10 — 100. We observed that non- 
axisymmetric motions in the disc lead to the formation of a 
radial component of the magnetic field, its stretching by the 
differential rotation in the disc, and subsequent MRI- driven 
turbulence in the disc. Interaction between the planet and 
turbulent cells leads to stochastic migration, similar to that 
observed in 3D discs by, e.g., [Nelson & Papaloizou] (|2004). 
We investigated the transition from the laminar to turbulent 
disc that often starts near vicinity of the planet, where non- 
axisymmetric density waves are excited by the planet. The 
turbulence starts more rapidly when the planet is more mas¬ 
sive. 

(v) The torques acting on the planet are larger in turbulent 
discs than in the laminar discs, leading to more rapid migra¬ 
tion. However, the direction of the migration is also stochas¬ 
tic, and outward migration is frequently observed. 

(vi) To understand how planets interact with individual 
turbulent cells, we investigated the propagation of a planet 
through density waves in the disc. We observed that a planet 
experiences a strong negative torque when it is located inside 
the wave or at the outer edge of the wave (where the den¬ 
sity increases towards the star). However, when the planet is 
located at the inner edge of the wave (where the density de¬ 
creases towards the star), then it experiences a positive coro¬ 
tation torque, and it migrates away from the star. We conclude 


that the stochastic motion of the planet is connected with the 
alternating action of these positive and negative torques. 
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